##### Functions to solve the firm first-order conditions

	function foc!(foc_values,base_premium)
		#CSV.write(iteration_file,DataFrame(cat(dims=2,plan_year_names,base_premium)));
		foc_moment_values = update_foc_values(base_premium,pmt,hdata,X_RS,X_CLMS,foc_residuals,nora_flag,mand_repeal_flag,voucher_flag,pc_flag,no_inertia_flag,no_uninsured_flag,ef_approach_flag);																							
		foc_values[1:length(plan_year_names)] = foc_moment_values; 
	end
	
	# Update Previous Choice Variable for Forward Simulation

	function update_prev_choice(base_premium,X_prev,households_prev,mand_repeal_flag,voucher_flag,no_uninsured_flag)
	
		# Initialize key variables
		X_demand = copy(X_prev);
		subsidy = voucher;
		new_unsubs_premiums = base_premium[household_pt_indices_prev] .* convert(Array,data_prev[uninsured_plans_prev .== 0,:age_geog_factor]);
		
		# Create data frame to hold objects we will group by household
		silverdf = DataFrame(household_number = household_numbers_prev[silver_indicators_prev],premium = new_unsubs_premiums[silver_indicators_prev],silver_indices = silver_indicators_prev);
		
		# Update Subsidy
		
			# Determine SLC plans by household
			num_household_silver_plans = combine(groupby(silverdf,:household_number), :household_number => length)[!,:household_number_length];
			pop!(num_household_silver_plans)
			num_household_silver_plans = cat(dims=1,0,num_household_silver_plans);
			starting_indices = cumsum(num_household_silver_plans); # starting position for each household
			household_slc_indices = combine(groupby(silverdf,:household_number), :premium => determine_SLC_plan)[!,:premium_determine_SLC_plan];
			slc_indices = silverdf[starting_indices .+ household_slc_indices,:silver_indices]; # SLC indices in hdata
			slc_plans = zeros(length(household_numbers_prev));
			slc_plans[slc_indices] = (households_prev[!,:subsidized_members] .> 0); # SLC index gets 1 if house is subsidized, 0 otherwise
		
			# Determine new subsidy (unless voucher)
			if !voucher_flag
				subsidy = max.(0,new_unsubs_premiums[slc_indices] - households_prev[!,:SLC_contribution]) .* 
					(households_prev[!,:subsidized_members] .> 0); # each household's subsidy (0 if unsubsidized)
			end
			subsidy = subsidy[household_numbers_prev];
			
		# Update Premium Variables
		X_demand[uninsured_plans_prev .== 0,findall(all_covariates .== "Premium")[1]] = 
				(max.(0,new_unsubs_premiums .- subsidy) .- 
					(1 - mand_repeal_flag) .* convert(Array,data_prev[uninsured_plans_prev .== 0,:penalty])) ./
				households_prev[household_numbers_prev,:enrollees];			
		X_demand[uninsured_plans_prev .== 0,premium_indices] = 
				X_demand[uninsured_plans_prev .== 0,premium_indices] .*  
				X_demand[uninsured_plans_prev .== 0,findall(all_covariates .== "Premium")[1]];
		
		if random_parameters
			probabilities = zeros(length(household_numbers_prev));
			ex_probs = zeros(length(household_numbers_prev));
			for n in 1:ns
				X_demand[:,random_parameter_indices] = X_prev[:,random_parameter_indices] .* V_prev[:,:,n]; 
				
				expdf = DataFrame(household_number = household_numbers_prev,
						expvarall = exp.((X_demand[uninsured_plans_prev .== 0,:] * (min_beta .* (1 .- previous_choice_indicator .* no_inertia_flag)))/lambda));
				if except_network_flag
					expdf[!,:expvarall] = exp.((X_demand[uninsured_plans_prev .== 0,:] * (min_beta .* (1 .- previous_choice_indicator_net .* except_network_flag)))/lambda);
				end
				expsum = max.(combine(groupby(expdf,:household_number), :expvarall => sum)[!,:expvarall_sum],0.0000000000001);
				
				if inattention_flag
					
					entrants = findall(households_prev[!,:previous_plan_number] .== "NA");
					search_probabilities = 1 ./ (1 .+ exp.(Xsearch_prev * beta_search));
					search_probabilities[entrants] .= 1; 
					
					uninsured_probabilities = 1  ./ (1 .+ expsum);
					insured_probabilities = search_probabilities[household_numbers_prev] .* expdf[!,:expvarall] ./ (1 .+ expsum[household_numbers_prev]) .+ 
						(1 .- search_probabilities[household_numbers_prev]) .* data_prev[uninsured_plans_prev .== 0,:previous_choice];
				else
					uninsured_probabilities = 1  ./ (1 .+ (expsum).^(lambda));
					insured_probabilities = expdf[!,:expvarall] .* (expsum[household_numbers_prev]).^(lambda-1) ./ (1 .+ (expsum[household_numbers_prev]).^(lambda));
				end
				
				probabilities += insured_probabilities ./ ns;
				ex_probs += (insured_probabilities ./ (1 .- uninsured_probabilities[household_numbers_prev])) ./ ns;
			end
		else
		
			expdf = DataFrame(household_number = household_numbers_prev,
						expvarall = exp.((X_demand[uninsured_plans_prev .== 0,:] * min_beta)/lambda));
			if except_network_flag
				expdf[!,:expvarall] = exp.((X_demand[uninsured_plans_prev .== 0,:] * (min_beta .* (1 .- previous_choice_indicator_net .* except_network_flag)))/lambda);
			end
			
			# Update choice probabilities
			expsum = combine(groupby(expdf,:household_number), :expvarall => sum)[!,:expvarall_sum];
			uninsured_probabilities = 1  ./ (1 .+ (expsum).^(lambda));
			probabilities = expdf[!,:expvarall] .* (expsum[household_numbers_prev]).^(lambda-1) ./ (1 .+ (expsum[household_numbers_prev]).^(lambda));
			ex_probs = probabilities ./ (1 .- uninsured_probabilities[household_numbers_prev]);
			
		end
		
		if no_uninsured_flag
			probabilities = ex_probs;
		end
		
		# Simulate new previous choices
		new_household_previous_choices = households_prev[!,:channel];
		full_probabilities = zeros(size(X_prev)[1]);
		full_probabilities[uninsured_plans_prev .== 0] = probabilities;
		full_probabilities[uninsured_plans_prev .== 1] = uninsured_probabilities;
		data_prev[!,:prob] = full_probabilities;
		
		#for i in 1:size(households_prev)[1]
		#	ids = findall(household_numbers_data_prev .== i);
		#	new_household_previous_choices[i] = sample(data_prev[ids,:plan_name],Weights(full_probabilities[ids]));
		#end
		
		new_household_previous_choices = combine(groupby(data_prev, :household_number), [:plan_name, :prob] => ((x, y) -> sample(x, Weights(y))) => :simulated_previous_choice);
		new_household_previous_choices = vec(new_household_previous_choices[!,2]);
		
		new_previous_choices = zeros(size(X)[1]);
		new_previous_choices[findall(in(string.(new_household_previous_choices,households_prev[!,:household_id])),
			string.(data[!,:plan_name],data[!,:household_id]))] .= 1;
		new_previous_choices[uninsured_plans .== 1] .= 0;
		
		if inattention_flag
			premiums_to_import = zeros(size(X_demand)[1]);
			premiums_to_import[uninsured_plans_prev .== 0] = new_unsubs_premiums;
			premiums_to_import = cat(dims=1,premiums_to_import,0);
			previous_choice_premiums = premiums_to_import[prev_prem_indices];
		end
		
		if inattention_flag
			return(cat(dims=2,new_previous_choices,previous_choice_premiums))
		else
			return(new_previous_choices)
		end
	end

	# Updates the previous choice variable and all interactions

	function update_previous_choice_variables(forward_simulated_prev_choices,X,households)
	
		X_demand = copy(X);
	
		X_demand[:,findall(all_covariates .== "previous_choice")[1]] = forward_simulated_prev_choices;
		
		X_demand[:,findall(all_covariates .== "previous_choice_250to400")[1]] = forward_simulated_prev_choices .* 
			((households[household_numbers_data,:FPL] .> 2.50) .& (households[household_numbers_data,:FPL] .<= 4)) ; 
		X_demand[:,findall(all_covariates .== "previous_choice_gt400")[1]] = forward_simulated_prev_choices .* 
			(households[household_numbers_data,:FPL] .> 4) ; 
		X_demand[:,findall(all_covariates .== "previous_choice_0to34")[1]] = 
			forward_simulated_prev_choices .* (households[household_numbers_data,:perc_0to17] .+ households[household_numbers_data,:perc_18to25] .+ 
				households[household_numbers_data,:perc_26to34]);
		X_demand[:,findall(all_covariates .== "previous_choice_35to54")[1]] = 
			forward_simulated_prev_choices .* (households[household_numbers_data,:perc_35to44] + households[household_numbers_data,:perc_45to54]);
		X_demand[:,findall(all_covariates .== "previous_choice_male")[1]] = 
			forward_simulated_prev_choices .* (households[household_numbers_data,:perc_male]);
		X_demand[:,findall(all_covariates .== "previous_choice_family")[1]] = 
			forward_simulated_prev_choices .* (households[household_numbers_data,:household_size] .> 1);
		X_demand[:,findall(all_covariates .== "previous_choice_Asian")[1]] = 
			forward_simulated_prev_choices .* (households[household_numbers_data,:perc_asian]);
		X_demand[:,findall(all_covariates .== "previous_choice_Black")[1]] = 
			forward_simulated_prev_choices .* (households[household_numbers_data,:perc_black]);
		X_demand[:,findall(all_covariates .== "previous_choice_Hispanic")[1]] = 
			forward_simulated_prev_choices .* (households[household_numbers_data,:perc_hispanic]);
		X_demand[:,findall(all_covariates .== "previous_choice_Other")[1]] = 
			forward_simulated_prev_choices .* (households[household_numbers_data,:perc_other]);
			
		X_demand[:,findall(all_covariates .== "previous_choice_250to400_0to34")[1]] = forward_simulated_prev_choices .* 
			((households[household_numbers_data,:FPL] .> 2.50) .& (households[household_numbers_data,:FPL] .<= 4)) .*
			(households[household_numbers_data,:perc_0to17] .+ households[household_numbers_data,:perc_18to25] .+ households[household_numbers_data,:perc_26to34]); 
		X_demand[:,findall(all_covariates .== "previous_choice_gt400_0to34")[1]] = forward_simulated_prev_choices .* 
			(households[household_numbers_data,:FPL] .> 4) .* (households[household_numbers_data,:perc_0to17] .+ households[household_numbers_data,:perc_18to25] .+ households[household_numbers_data,:perc_26to34]);  
		X_demand[:,findall(all_covariates .== "previous_choice_250to400_35to54")[1]] = forward_simulated_prev_choices .* 
			((households[household_numbers_data,:FPL] .> 2.50) .& (households[household_numbers_data,:FPL] .<= 4)) .* (households[household_numbers_data,:perc_35to44] .+ households[household_numbers_data,:perc_45to54]); 
		X_demand[:,findall(all_covariates .== "previous_choice_gt400_35to54")[1]] = forward_simulated_prev_choices .* 
			(households[household_numbers_data,:FPL] .> 4) .* (households[household_numbers_data,:perc_35to44] .+ households[household_numbers_data,:perc_45to54]);  
				
			
		#X_demand[:,findall(all_covariates .== "previous_choice_Anthem")[1]] = 
		#	forward_simulated_prev_choices .* X_demand[:,findall(all_covariates .== "Anthem")[1]];
		#X_demand[:,findall(all_covariates .== "previous_choice_Blue_Shield")[1]] = 
		#	forward_simulated_prev_choices .* X_demand[:,findall(all_covariates .== "Blue_Shield")[1]];
		#X_demand[:,findall(all_covariates .== "previous_choice_Kaiser")[1]] = 
		#	forward_simulated_prev_choices .* X_demand[:,findall(all_covariates .== "Kaiser")[1]];
		#X_demand[:,findall(all_covariates .== "previous_choice_Health_Net")[1]] = 
		#	forward_simulated_prev_choices .* X_demand[:,findall(all_covariates .== "Health_Net")[1]];			
		#X_demand[:,findall(all_covariates .== "previous_choice_HMO")[1]] = 
		#	forward_simulated_prev_choices .* X_demand[:,findall(all_covariates .== "HMO")[1]];
		#X_demand[:,findall(all_covariates .== "previous_choice_AV")[1]] = 
		#	forward_simulated_prev_choices .* X_demand[:,findall(all_covariates .== "AV")[1]];
		#if add_silver
		#	X_demand[:,findall(all_covariates .== "previous_choice_silver")[1]] = 
		#		forward_simulated_prev_choices .* X_demand[:,findall(all_covariates .== "Silver")[1]];
		#end
		if network_flag 
			X_demand[:,findall(all_covariates .== "previous_choice_breadth")[1]] = 
				forward_simulated_prev_choices .* convert(Array,data[!,:network_breadth]);
			X_demand[:,findall(all_covariates .== "previous_choice_inclus")[1]] = 
				forward_simulated_prev_choices .* convert(Array,data[!,:network_inclus]);
		end
	
		return(X_demand)
	end

	# Update All Variables for Counterfactual

	function update_foc_values(base_premium,pmt,hdata,X_RS,X_CLMS,foc_residuals,nora_flag,mand_repeal_flag,voucher_flag,pc_flag,no_inertia_flag,no_uninsured_flag,ef_approach_flag)
		
		# Initialize key variables
		X_demand = copy(X);
		Xsearch_df = copy(Xsearch);
		subsidy = voucher;
		age_factor = hdata[:,findall(hdata_fields .== "age_factor")[1]];
		rating_factor = hdata[:,findall(hdata_fields .== "rating_factor")[1]];
		new_unsubs_premiums = base_premium[household_pt_indices] .* rating_factor;
		
		# Create data frame to hold objects we will group by household
		silverdf = DataFrame(household_number = household_numbers[silver_indicators],premium = new_unsubs_premiums[silver_indicators],silver_indices = silver_indicators);
		
		# Update Subsidy
		
			# Determine SLC plans by household
			num_household_silver_plans = combine(groupby(silverdf,:household_number), :household_number => length)[!,:household_number_length];
			pop!(num_household_silver_plans)
			num_household_silver_plans = cat(dims=1,0,num_household_silver_plans);
			starting_indices = cumsum(num_household_silver_plans); # starting position for each household
			household_slc_indices = combine(groupby(silverdf,:household_number), :premium => determine_SLC_plan)[!,:premium_determine_SLC_plan];
			slc_indices = silverdf[starting_indices .+ household_slc_indices,:silver_indices]; # SLC indices in hdata
			slc_plans = zeros(size(hdata)[1]);
			slc_plans[slc_indices] = (households[!,:subsidized_members] .> 0); # SLC index gets 1 if house is subsidized, 0 otherwise
		
			# Determine new subsidy (unless voucher)
			if !voucher_flag
				subsidy = max.(0,new_unsubs_premiums[slc_indices] - households[!,:SLC_contribution]) .* 
					(households[!,:subsidized_members] .> 0); # each household's subsidy (0 if unsubsidized)
				subsidy = subsidy[household_numbers];
			end
		
		# Update Premium Variables
		X_demand[uninsured_plans .== 0,findall(all_covariates .== "Premium")[1]] = 
				(max.(0,new_unsubs_premiums .- subsidy) .- 
					(1 - mand_repeal_flag) .* hdata[:,findall(hdata_fields .== "penalty")[1]]) ./
				hdata[:,findall(hdata_fields .== "members")[1]];			
		X_demand[uninsured_plans .== 0,premium_indices] = 
				X_demand[uninsured_plans .== 0,premium_indices] .*  
				X_demand[uninsured_plans .== 0,findall(all_covariates .== "Premium")[1]];
		
		if inattention_flag & (year > 2014) & search_flag
		
			# Growth rates
			growth_rates = new_unsubs_premiums ./ forward_simulated_prev_premiums[uninsured_plans .== 0];
			fsim_prev_choices = forward_simulated_prev_choices[uninsured_plans .== 0];
			growthdf = DataFrame(household_number = household_numbers[offered_indices],growth_rate = growth_rates[offered_indices],prev_choice=fsim_prev_choices[offered_indices]);
			median_growth_rates = combine(groupby(growthdf,:household_number), :growth_rate => median);
			
			household_median_growth_rates = zeros(size(households)[1]);
			household_median_growth_rates[median_growth_rates[!,:household_number]] = median_growth_rates[!,:growth_rate_median];
			
			growthdf_for_choices = growthdf[findall(growthdf[!,:prev_choice] .== 1),:];
			growthdf_for_choices[!,:median_growth_rate] = household_median_growth_rates[growthdf_for_choices[!,:household_number]];
			
			premium_shock_households = growthdf_for_choices[findall(growthdf_for_choices[!,:growth_rate] .> growthdf_for_choices[!,:median_growth_rate]),:household_number];
			Xsearch_df[:,findall(search_covariates .== "premium_shock")[1]] .= 0;
			Xsearch_df[premium_shock_households,findall(search_covariates .== "premium_shock")[1]] .= 1;
		
		end
		
		if random_parameters
			#X_demand[:,random_parameter_indices] = X_demand[:,random_parameter_indices] .* V[:,:,1]; 
			probabilities = zeros(size(hdata)[1]);
			ex_probs = zeros(size(hdata)[1]);
			#uni_probs = zeros(size(households)[1]);
			for n in 1:ns
				X_demand[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
				
				expdf = DataFrame(household_number = household_numbers,
						expvarall = exp.((X_demand[uninsured_plans .== 0,:] * (min_beta .* (1 .- previous_choice_indicator .* no_inertia_flag)))/lambda));
				if except_network_flag
					expdf[!,:expvarall] = exp.((X_demand[uninsured_plans .== 0,:] * (min_beta .* (1 .- previous_choice_indicator_net .* except_network_flag)))/lambda);
				end
				expsum = max.(combine(groupby(expdf,:household_number), :expvarall => sum)[!,:expvarall_sum],0.0000000000001);
				
				if inattention_flag & search_flag & (year > 2014)
					entrants = findall(households[!,:previous_plan_number] .== "NA");
					search_probabilities = 1 ./ (1 .+ exp.(Xsearch_df * beta_search));
					search_probabilities[entrants] .= 1; 
					
					uninsured_probabilities = 1  ./ (1 .+ expsum);
					insured_probabilities = search_probabilities[household_numbers] .* expdf[!,:expvarall] ./ (1 .+ expsum[household_numbers]) .+ 
						(1 .- search_probabilities[household_numbers]) .* data[uninsured_plans .== 0,:previous_choice];
				else
					uninsured_probabilities = 1  ./ (1 .+ (expsum).^(lambda));
					insured_probabilities = expdf[!,:expvarall] .* (expsum[household_numbers]).^(lambda-1) ./ (1 .+ (expsum[household_numbers]).^(lambda));
				end
				
				probabilities += insured_probabilities ./ ns;
				ex_probs += (insured_probabilities ./ (1 .- uninsured_probabilities[household_numbers])) ./ ns;
				#uni_probs += uninsured_probabilities ./ ns;
			end
		else
		
			expdf = DataFrame(household_number = household_numbers,
						expvarall = exp.((X_demand[uninsured_plans .== 0,:] * (min_beta .* (1 .- previous_choice_indicator .* no_inertia_flag)))/lambda));
			if except_network_flag
				expdf[!,:expvarall] = exp.((X_demand[uninsured_plans .== 0,:] * (min_beta .* (1 .- previous_choice_indicator_net .* except_network_flag)))/lambda);
			end
			
			# Update choice probabilities
			expsum = combine(groupby(expdf,:household_number), :expvarall => sum)[!,:expvarall_sum];
			uninsured_probabilities = 1  ./ (1 .+ (expsum).^(lambda));
			probabilities = expdf[!,:expvarall] .* (expsum[household_numbers]).^(lambda-1) ./ (1 .+ (expsum[household_numbers]).^(lambda));
			ex_probs = probabilities ./ (1 .- uninsured_probabilities[household_numbers]);
			
		end
		
		if no_uninsured_flag
			probabilities = ex_probs;
		end
		
		# Calculate weights
		probability_weights = probabilities .* hdata[:,findall(hdata_fields .== "weight")[1]];
		pricing_weights = probabilities .* age_factor;
		age_weights = probability_weights .* (age_factor ./ hdata[:,findall(hdata_fields .== "members")[1]]);
		rf_pricing_weights = probabilities .* rating_factor;
		if old_flag
			weight_object = DataFrame(household_number = household_numbers, plan_name = data_pmt_names,
				probability_weights = probability_weights,pricing_weights = pricing_weights,rf_pricing_weights = rf_pricing_weights,age_weights = age_weights,
				rs_share_weights_18to54 = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_18to54")[1]],
				rs_share_weights_Hispanic = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_Hispanic")[1]]);
		elseif agegender_flag
			weight_object = DataFrame(household_number = household_numbers, plan_name = data_pmt_names,
				probability_weights = probability_weights,pricing_weights = pricing_weights,rf_pricing_weights = rf_pricing_weights,age_weights = age_weights,
				rs_share_weights_0to34 = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_0to34")[1]],
				rs_share_weights_male = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_male")[1]]);
		elseif hisp_flag
			weight_object = DataFrame(household_number = household_numbers, plan_name = data_pmt_names,
				probability_weights = probability_weights,pricing_weights = pricing_weights,rf_pricing_weights = rf_pricing_weights,age_weights = age_weights,
				rs_share_weights_0to34 = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_0to34")[1]],
				rs_share_weights_male = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_male")[1]],
				rs_share_weights_family = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_family")[1]],
				rs_share_weights_Hispanic = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_Hispanic")[1]]);
		elseif int_rs_variables
			weight_object = DataFrame(household_number = household_numbers, plan_name = data_pmt_names,
				probability_weights = probability_weights,pricing_weights = pricing_weights,rf_pricing_weights = rf_pricing_weights,age_weights = age_weights,
				rs_share_weights_0to34 = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_0to34")[1]],
				rs_share_weights_male = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_male")[1]],
				rs_share_weights_family = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_family")[1]],
				rs_share_weights_minority = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_minority")[1]]);
		elseif all_rs_variables
			weight_object = DataFrame(household_number = household_numbers, plan_name = data_pmt_names,
				probability_weights = probability_weights,pricing_weights = pricing_weights,rf_pricing_weights = rf_pricing_weights,age_weights = age_weights,
				rs_share_weights_0to34 = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_0to34")[1]],
				rs_share_weights_35to54 = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_35to54")[1]],
				rs_share_weights_250to400 = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_250to400")[1]],
				rs_share_weights_gt400 = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_gt400")[1]],
				rs_share_weights_male = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_male")[1]],
				rs_share_weights_family = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_family")[1]],
				rs_share_weights_Asian = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_Asian")[1]],
				rs_share_weights_Black = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_Black")[1]],
				rs_share_weights_Hispanic = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_Hispanic")[1]],
				rs_share_weights_Other = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_Other_Race")[1]]);
		end
		
		# Update plan-level variables
		plan_pricing_weights = sort!(combine(groupby(weight_object,:plan_name), :pricing_weights => sum),[:plan_name])[!,:pricing_weights_sum];
		plan_demands = sort!(combine(groupby(weight_object,:plan_name), :probability_weights => sum),[:plan_name])[!,:probability_weights_sum];
		plan_total_premiums = base_premium[year_plan_indices] .* 
			sort!(combine(groupby(weight_object,:plan_name), :rf_pricing_weights => sum),[:plan_name])[!,:rf_pricing_weights_sum];
		plan_average_age_factors = 
			sort!(combine(groupby(weight_object,:plan_name), :age_weights => sum),[:plan_name])[!,:age_weights_sum] ./ plan_demands;	
		
		# Update risk score variables
		X_RS_new = copy(X_RS);
		if old_flag
			X_RS_new[:,findall(risk_score_covariates .== "share_18to54")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_18to54 => sum),[:plan_name])[!,:rs_share_weights_18to54_sum] ./ plan_demands;
			X_RS_new[:,findall(risk_score_covariates .== "share_Hispanic")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_Hispanic => sum),[:plan_name])[!,:rs_share_weights_Hispanic_sum] ./ plan_demands;
		elseif agegender_flag
			X_RS_new[:,findall(risk_score_covariates .== "share_0to34")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_0to34 => sum),[:plan_name])[!,:rs_share_weights_0to34_sum] ./ plan_demands;
			X_RS_new[:,findall(risk_score_covariates .== "share_male")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_male => sum),[:plan_name])[!,:rs_share_weights_male_sum] ./ plan_demands;
		elseif hisp_flag
			X_RS_new[:,findall(risk_score_covariates .== "share_0to34")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_0to34 => sum),[:plan_name])[!,:rs_share_weights_0to34_sum] ./ plan_demands;
			X_RS_new[:,findall(risk_score_covariates .== "share_male")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_male => sum),[:plan_name])[!,:rs_share_weights_male_sum] ./ plan_demands;
			X_RS_new[:,findall(risk_score_covariates .== "share_family")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_family => sum),[:plan_name])[!,:rs_share_weights_family_sum] ./ plan_demands;
			X_RS_new[:,findall(risk_score_covariates .== "share_Hispanic")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_Hispanic => sum),[:plan_name])[!,:rs_share_weights_Hispanic_sum] ./ plan_demands;
		elseif int_rs_variables
			X_RS_new[:,findall(risk_score_covariates .== "share_0to34")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_0to34 => sum),[:plan_name])[!,:rs_share_weights_0to34_sum] ./ plan_demands;
			X_RS_new[:,findall(risk_score_covariates .== "share_male")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_male => sum),[:plan_name])[!,:rs_share_weights_male_sum] ./ plan_demands;
			X_RS_new[:,findall(risk_score_covariates .== "share_family")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_family => sum),[:plan_name])[!,:rs_share_weights_family_sum] ./ plan_demands;
			X_RS_new[:,findall(risk_score_covariates .== "share_minority")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_minority => sum),[:plan_name])[!,:rs_share_weights_minority_sum] ./ plan_demands;
		elseif all_rs_variables
			X_RS_new[:,findall(risk_score_covariates .== "share_0to34")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_0to34 => sum),[:plan_name])[!,:rs_share_weights_0to34_sum] ./ plan_demands;
			X_RS_new[:,findall(risk_score_covariates .== "share_35to54")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_35to54 => sum),[:plan_name])[!,:rs_share_weights_35to54_sum] ./ plan_demands;
			X_RS_new[:,findall(risk_score_covariates .== "share_250to400")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_250to400 => sum),[:plan_name])[!,:rs_share_weights_250to400_sum] ./ plan_demands;
			X_RS_new[:,findall(risk_score_covariates .== "share_gt400")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_gt400 => sum),[:plan_name])[!,:rs_share_weights_gt400_sum] ./ plan_demands;
			X_RS_new[:,findall(risk_score_covariates .== "share_male")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_male => sum),[:plan_name])[!,:rs_share_weights_male_sum] ./ plan_demands;
			X_RS_new[:,findall(risk_score_covariates .== "share_family")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_family => sum),[:plan_name])[!,:rs_share_weights_family_sum] ./ plan_demands;
			X_RS_new[:,findall(risk_score_covariates .== "share_Asian")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_Asian => sum),[:plan_name])[!,:rs_share_weights_Asian_sum] ./ plan_demands;
			X_RS_new[:,findall(risk_score_covariates .== "share_Black")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_Black => sum),[:plan_name])[!,:rs_share_weights_Black_sum] ./ plan_demands;
			X_RS_new[:,findall(risk_score_covariates .== "share_Hispanic")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_Hispanic => sum),[:plan_name])[!,:rs_share_weights_Hispanic_sum] ./ plan_demands;
			X_RS_new[:,findall(risk_score_covariates .== "share_Other_Race")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_Other => sum),[:plan_name])[!,:rs_share_weights_Other_sum] ./ plan_demands;
		end
		
		plan_risk_scores = exp.(X_RS_new * theta_RS);
		household_risk_scores = plan_risk_scores[household_pmt_indices];	
			
		# Calculate cost factors
		rs_cost_factor = plans_pmt[!,:IDF] .* plans_pmt[!,:GCF] .* plan_risk_scores;
		ra_cost_factor = plans_pmt[!,:AV] .* plans_pmt[!,:IDF] .* plans_pmt[!,:GCF] .* plan_average_age_factors;
		household_ra_cost_factor = (age_factor ./ hdata[:,findall(hdata_fields .== "members")[1]]) .*
			hdata[:,findall(hdata_fields .== "ra_factor")[1]] .* hdata[:,findall(hdata_fields .== "GCF")[1]];		
			
		# Update household-level variables that depend on household jacobian
		jacobian = zeros(length(plan_pmt_names),length(plan_pmt_names));
		rs_partial_matrix = zeros(length(plan_pmt_names),length(plan_pmt_names));
		price_partial_firmprob_partial_price_matrix = zeros(length(insurers),length(plan_pmt_names));
		partial_rs_demographic_own_demand_partial_price = zeros(length(plan_pmt_names));
		partial_rs_demographic_all_demand_partial_price = zeros(length(plan_pmt_names));
		partial_ra_demographic_own_demand_partial_price = zeros(length(plan_pmt_names));
		partial_ra_demographic_all_demand_partial_price = zeros(length(plan_pmt_names));
			
		# Base vectors to form variables that depend on household-level jacobian
		cross_vector1 = (lambda-1)/lambda * ex_probs .- probabilities;
		cross_vector2 = probabilities .* hdata[:,findall(hdata_fields .== "alpha")[1]]
		own_vector = (1/lambda .+ (lambda-1)/lambda * ex_probs .- probabilities) .* probabilities .* hdata[:,findall(hdata_fields .== "alpha")[1]];
		household_premiums = rating_factor .* base_premium[household_pt_indices];
		
		# NOTE: You can make this loop more modular, but it slows down the code
		
		for m in 1:size(unique_choice_sets)[1]
			
			# Indices
		
				# Get data indices
				choice_set_indices = findall(data[uninsured_plans .== 0,:choice_set_id] .== m);
				subsidized_indices = intersect(choice_set_indices,findall(households[household_numbers,:subsidized_members] .> 0));
				unsubsidized_indices = intersect(choice_set_indices,findall(households[household_numbers,:subsidized_members] .== 0));
					
				# Get household indices	
				household_indices = unique(household_numbers[choice_set_indices]);
				subsidized_household_indices = unique(household_numbers[subsidized_indices]);
				unsubsidized_household_indices = unique(household_numbers[unsubsidized_indices]); 
				
				# Get plan indices
				jacobian_indices = unique(household_pmt_indices[choice_set_indices]);
				if length(subsidized_indices) > 0
					choice_slc_jacobian_index = findall(slc_plans[subsidized_indices] .> 0)[1];
					choice_nonslc_jacobian_indices = setdiff(1:length(jacobian_indices),choice_slc_jacobian_index);
				end
				firm_indices = findall(vec(sum(dims=2,firm_plan_matrix[:,jacobian_indices])) .> 0);	
				
				num_households_subs = length(subsidized_household_indices);
				num_households_unsubs = length(unsubsidized_household_indices);
				num_plans = length(jacobian_indices);
			
			# Matrices we will reuse many times to calculate all jacobian variants
			if num_households_unsubs > 0
				cross_matrix1_unsubs = reshape(cross_vector1[unsubsidized_indices],num_plans,num_households_unsubs);
				cross_matrix2_unsubs = reshape(cross_vector2[unsubsidized_indices],num_plans,num_households_unsubs);
				own_matrix_unsubs = reshape(own_vector[unsubsidized_indices],num_plans,num_households_unsubs);
				household_premiums_unsubs = reshape(household_premiums[unsubsidized_indices],num_plans,num_households_unsubs);
			end
			if num_households_subs > 0
				cross_matrix1_subs = reshape(cross_vector1[subsidized_indices],num_plans,num_households_subs);
				cross_matrix2_subs = reshape(cross_vector2[subsidized_indices],num_plans,num_households_subs);
				own_matrix_subs = reshape(own_vector[subsidized_indices],num_plans,num_households_subs);
				household_premiums_subs = reshape(household_premiums[subsidized_indices],num_plans,num_households_subs);
			end	
			
			# Get price jacobian
			
				if num_households_unsubs > 0
					multiplier_unsubs = reshape(rating_factor[unsubsidized_indices] .* household_weights[unsubsidized_indices],num_plans,num_households_unsubs);
					unsubsidized_portion = cross_matrix1_unsubs * transpose(cross_matrix2_unsubs .* multiplier_unsubs);
					unsubsidized_portion[diagind(unsubsidized_portion)] = sum(dims=2,own_matrix_unsubs .* multiplier_unsubs);
					jacobian[jacobian_indices,jacobian_indices] += unsubsidized_portion;
				end
			
				if num_households_subs > 0
					multiplier_subs = reshape(rating_factor[subsidized_indices] .* household_weights[subsidized_indices],num_plans,num_households_subs);
					subsidized_portion = cross_matrix1_subs * transpose(cross_matrix2_subs .* multiplier_subs);
					subsidized_portion[diagind(subsidized_portion)] = sum(dims=2,own_matrix_subs .* multiplier_subs);
					
					# Apply adjustment for endogenous subsidy (only affects SLC column of subsidized portion)
							# want to calculate partial q_l/partial p_SLC for all l (including SLC)
							# negative sum of all partial q_l/partial p_l', except l' = SLC 
					
					if !voucher_flag
						multiplier_slc = reshape(repeat(rating_factor[intersect(subsidized_indices,slc_indices)],
							inner=length(jacobian_indices)) .* household_weights[subsidized_indices],num_plans,num_households_subs);
						slc_portion = cross_matrix1_subs * transpose(cross_matrix2_subs .* multiplier_slc);
						slc_portion[diagind(slc_portion)] = sum(dims=2,own_matrix_subs .* multiplier_slc);
						subsidized_portion[:,choice_slc_jacobian_index] = -vec(sum(dims=2,slc_portion[:,choice_nonslc_jacobian_indices]));
					end
					
					jacobian[jacobian_indices,jacobian_indices] += subsidized_portion;
				end
			
			# Populate the risk score partial matrix (input for the claim slope matrix)
					# Note this is a portion; we'll subtract the second portion later
			
				if num_households_unsubs > 0
					rs_multiplier_unsubs = reshape((rs_cost_factor[jacobian_indices] ./ plan_demands[jacobian_indices]) * 
						transpose(hh_share_matrix[unsubsidized_household_indices,:] * theta_RS_share_variables),num_plans,num_households_unsubs);
					unsubsidized_portion = (cross_matrix1_unsubs .* rs_multiplier_unsubs) * transpose(cross_matrix2_unsubs .* multiplier_unsubs);
					unsubsidized_portion[diagind(unsubsidized_portion)] = sum(dims=2,own_matrix_unsubs .* multiplier_unsubs .* rs_multiplier_unsubs);
					rs_partial_matrix[jacobian_indices,jacobian_indices] += unsubsidized_portion;
				end
			
				if num_households_subs > 0
					rs_multiplier_subs = reshape((rs_cost_factor[jacobian_indices] ./ plan_demands[jacobian_indices]) * 
						transpose(hh_share_matrix[subsidized_household_indices,:] * theta_RS_share_variables),num_plans,num_households_subs);
					subsidized_portion = (cross_matrix1_subs .* rs_multiplier_subs) * transpose(cross_matrix2_subs .* multiplier_subs);
					subsidized_portion[diagind(subsidized_portion)] = sum(dims=2,own_matrix_subs .* multiplier_subs .* rs_multiplier_subs);
					
					if !voucher_flag					
						slc_portion = (cross_matrix1_subs .* rs_multiplier_subs) * transpose(cross_matrix2_subs .* multiplier_slc);
						slc_portion[diagind(slc_portion)] = sum(dims=2,own_matrix_subs .* multiplier_slc .* rs_multiplier_subs);
						subsidized_portion[:,choice_slc_jacobian_index] = -vec(sum(dims=2,slc_portion[:,choice_nonslc_jacobian_indices])); 	
					end
					
					rs_partial_matrix[jacobian_indices,jacobian_indices] += subsidized_portion;						
				end
			
			# Get partials that go into marginal RA transfer

				if num_households_unsubs > 0
					rs_multiplier_unsubs = reshape(rs_cost_factor[jacobian_indices] * 
						transpose(hh_share_matrix[unsubsidized_household_indices,:] * theta_RS_share_variables),num_plans,num_households_unsubs);
					unsubsidized_portion = (cross_matrix1_unsubs .* rs_multiplier_unsubs) * transpose(cross_matrix2_unsubs .* multiplier_unsubs);
					unsubsidized_portion[diagind(unsubsidized_portion)] = sum(dims=2,own_matrix_unsubs .* multiplier_unsubs .* rs_multiplier_unsubs);
				
					partial_rs_demographic_own_demand_partial_price[jacobian_indices] += vec(sum(dims=1,unsubsidized_portion .* own_plan_year_matrix[jacobian_indices,jacobian_indices]));
					partial_rs_demographic_all_demand_partial_price[jacobian_indices] += vec(sum(dims=1,unsubsidized_portion));
					
					ra_multiplier_unsubs = reshape(household_ra_cost_factor[unsubsidized_indices],num_plans,num_households_unsubs);;
					unsubsidized_portion = (cross_matrix1_unsubs .* ra_multiplier_unsubs) * transpose(cross_matrix2_unsubs .* multiplier_unsubs);
					unsubsidized_portion[diagind(unsubsidized_portion)] = sum(dims=2,own_matrix_unsubs .* multiplier_unsubs .* ra_multiplier_unsubs);
					
					partial_ra_demographic_own_demand_partial_price[jacobian_indices] += vec(sum(dims=1,unsubsidized_portion .* own_plan_year_matrix[jacobian_indices,jacobian_indices]));
					partial_ra_demographic_all_demand_partial_price[jacobian_indices] += vec(sum(dims=1,unsubsidized_portion));
				
				end

				if num_households_subs > 0
					rs_multiplier_subs = reshape(rs_cost_factor[jacobian_indices] * 
						transpose(hh_share_matrix[subsidized_household_indices,:] * theta_RS_share_variables),num_plans,num_households_subs);
					subsidized_portion = (cross_matrix1_subs .* rs_multiplier_subs) * transpose(cross_matrix2_subs .* multiplier_subs);
					subsidized_portion[diagind(subsidized_portion)] = sum(dims=2,own_matrix_subs .* multiplier_subs .* rs_multiplier_subs);
					
					if !voucher_flag
						slc_portion = (cross_matrix1_subs .* rs_multiplier_subs) * transpose(cross_matrix2_subs .* multiplier_slc);
						slc_portion[diagind(slc_portion)] = sum(dims=2,own_matrix_subs .* multiplier_slc .* rs_multiplier_subs);
						subsidized_portion[:,choice_slc_jacobian_index] = -vec(sum(dims=2,slc_portion[:,choice_nonslc_jacobian_indices]));
					end
					
					partial_rs_demographic_own_demand_partial_price[jacobian_indices] += vec(sum(dims=1,subsidized_portion .* own_plan_year_matrix[jacobian_indices,jacobian_indices]));
					partial_rs_demographic_all_demand_partial_price[jacobian_indices] += vec(sum(dims=1,subsidized_portion));

					ra_multiplier_subs = reshape(household_ra_cost_factor[subsidized_indices],num_plans,num_households_subs);;
					subsidized_portion = (cross_matrix1_subs .* ra_multiplier_subs) * transpose(cross_matrix2_subs .* multiplier_subs);
					subsidized_portion[diagind(subsidized_portion)] = sum(dims=2,own_matrix_subs .* multiplier_subs .* ra_multiplier_subs);

					if !voucher_flag
						slc_portion = (cross_matrix1_subs .* ra_multiplier_subs) * transpose(cross_matrix2_subs .* multiplier_slc);
						slc_portion[diagind(slc_portion)] = sum(dims=2,own_matrix_subs .* multiplier_slc .* ra_multiplier_subs);
						subsidized_portion[:,choice_slc_jacobian_index] = -vec(sum(dims=2,slc_portion[:,choice_nonslc_jacobian_indices]));
					end
					
					partial_ra_demographic_own_demand_partial_price[jacobian_indices] += vec(sum(dims=1,subsidized_portion .* own_plan_year_matrix[jacobian_indices,jacobian_indices]));
					partial_ra_demographic_all_demand_partial_price[jacobian_indices] += vec(sum(dims=1,subsidized_portion));
				end

			# Get partial_firmprob_partial_price (do not multiply by household's weight)
			
				if num_households_unsubs > 0
					multiplier_unsubs = reshape(rating_factor[unsubsidized_indices],num_plans,num_households_unsubs);
					unsubsidized_portion = (cross_matrix1_unsubs .* household_premiums_unsubs) * transpose(cross_matrix2_unsubs .* multiplier_unsubs);
					unsubsidized_portion[diagind(unsubsidized_portion)] = sum(dims=2,own_matrix_unsubs .* multiplier_unsubs .* household_premiums_unsubs);
					price_partial_firmprob_partial_price_matrix[firm_indices,jacobian_indices] += 
						firm_plan_matrix[firm_indices,jacobian_indices] * unsubsidized_portion;
				end

				if num_households_subs > 0
					multiplier_subs = reshape(rating_factor[subsidized_indices],num_plans,num_households_subs);
					subsidized_portion = (cross_matrix1_subs .* household_premiums_subs) * transpose(cross_matrix2_subs .* multiplier_subs);
					subsidized_portion[diagind(subsidized_portion)] = sum(dims=2,own_matrix_subs .* multiplier_subs .* household_premiums_subs);
					
					if !voucher_flag
						choice_slc_indices = intersect(subsidized_indices,slc_indices);
						multiplier_slc = reshape(repeat(rating_factor[choice_slc_indices],inner=length(jacobian_indices)),num_plans,num_households_subs);
						slc_portion = (cross_matrix1_subs .* household_premiums_subs) * transpose(cross_matrix2_subs .* multiplier_slc);
						slc_portion[diagind(slc_portion)] = sum(dims=2,own_matrix_subs .* multiplier_slc .* household_premiums_subs);
						subsidized_portion[:,choice_slc_jacobian_index] = -vec(sum(dims=2,slc_portion[:,choice_nonslc_jacobian_indices])); 
					end
					
					price_partial_firmprob_partial_price_matrix[firm_indices,jacobian_indices] += 
						firm_plan_matrix[firm_indices,jacobian_indices] * subsidized_portion;
				end
		end	
		
		# Calculate marginal revenue
			
			partial_firmrevenue_partial_price_matrix = price_partial_firmprob_partial_price_matrix .+ firm_plan_matrix .* transpose(plan_pricing_weights);
			own_firmrevenue_indices = findall(transpose(firm_plan_matrix) .!= 0);			
			partial_firmrevenue_partial_price = transpose(partial_firmrevenue_partial_price_matrix)[own_firmrevenue_indices];
			partial_marketrevenue_partial_price = vec(sum(dims=1,partial_firmrevenue_partial_price_matrix));
			
		# Calculate marginal transfer
		
			annual_total_premiums = transpose(plan_year_matrix) * pmt[:,findall(pmt_fields .== "total_premiums")[1]];
			annual_demand = transpose(plan_year_matrix) * plan_demands;
			annual_ra_weights = transpose(plan_year_matrix) * (plan_demands .* ra_cost_factor);
			annual_rs_weights = transpose(plan_year_matrix) * (plan_demands .* rs_cost_factor);
			
			insurer_rs_weights = transpose(own_plan_year_matrix) * (plan_demands .* rs_cost_factor);
			insurer_ra_weights = transpose(own_plan_year_matrix) * (plan_demands .* ra_cost_factor);	
			insurer_rs_shares = insurer_rs_weights ./ annual_rs_weights;
			insurer_ra_shares = insurer_ra_weights ./ annual_ra_weights;
			
			rs_partial_ownplans_partial_price = transpose(jacobian .* own_plan_year_matrix) * rs_cost_factor;
			rs_partial_allplans_partial_price = transpose(jacobian) * rs_cost_factor;
			ra_partial_ownplans_partial_price = transpose(jacobian .* own_plan_year_matrix) * ra_cost_factor;
			ra_partial_allplans_partial_price = transpose(jacobian) * ra_cost_factor;
			
			partial_rs_dem_share_own_plan_demand_partial_price = transpose(jacobian .* own_plan_year_matrix) * 
				(rs_cost_factor .* (X_RS[:,theta_RS_share_indices] * theta_RS_share_variables));
			partial_rs_dem_share_all_plan_demand_partial_price = transpose(jacobian) * 
				(rs_cost_factor .* (X_RS[:,theta_RS_share_indices] * theta_RS_share_variables));		
	
			partial_ra_dem_share_own_plan_demand_partial_price = (transpose(jacobian .* own_plan_year_matrix) * ra_cost_factor);
			partial_ra_dem_share_all_plan_demand_partial_price = (transpose(jacobian) * ra_cost_factor);	
	
			own_plan_firmshare_partial_rs_partial_price =
				partial_rs_demographic_own_demand_partial_price .- partial_rs_dem_share_own_plan_demand_partial_price;
			own_plan_firmshare_partial_ra_partial_price =
				partial_ra_demographic_own_demand_partial_price .- partial_ra_dem_share_own_plan_demand_partial_price;
			
			all_plan_firmshare_partial_rs_partial_price =			
				partial_rs_demographic_all_demand_partial_price .- partial_rs_dem_share_all_plan_demand_partial_price;
			all_plan_firmshare_partial_ra_partial_price =			
				partial_ra_demographic_all_demand_partial_price .- partial_ra_dem_share_all_plan_demand_partial_price;
			
			partial_firmriskshare_partial_price = (1 ./ annual_rs_weights) .*
				((own_plan_firmshare_partial_rs_partial_price .+ rs_partial_ownplans_partial_price) .-
					insurer_rs_shares .* (all_plan_firmshare_partial_rs_partial_price .+ rs_partial_allplans_partial_price));	
			partial_firmutilizationshare_partial_price = (1 ./ annual_ra_weights) .*
				((own_plan_firmshare_partial_ra_partial_price .+ ra_partial_ownplans_partial_price) .-
					insurer_ra_shares .* (all_plan_firmshare_partial_ra_partial_price .+ ra_partial_allplans_partial_price));	
			
			partial_firmtransfer_partial_price = nu_for_riskadj .*
				(annual_total_premiums .* (partial_firmriskshare_partial_price .- partial_firmutilizationshare_partial_price) .+
				partial_marketrevenue_partial_price .* (insurer_rs_shares .- insurer_ra_shares));
			
		# Calculate second part of risk score partial matrix
			
			rs_partial_matrix -= jacobian .*
				(rs_cost_factor .* (X_RS_new[:,theta_RS_share_indices] * theta_RS_share_variables) ./ plan_demands);
		
		# Calculate marginal claims
		
			X_CLMS_new = copy(X_CLMS);
			X_CLMS_new[:,findall(claims_moments_covariates .== "log_risk_score")[1]] = X_RS_new * theta_RS;
			average_claims = exp.(X_CLMS_new * theta_CLMS);
			
			ac_partial_ownplans_partial_price = transpose(jacobian .* own_plan_year_matrix) * average_claims;
			demand_partial_ac_partial_price = theta_CLMS[findall(claims_moments_covariates .== "log_risk_score")[1]] .*
				(transpose(rs_partial_matrix .* own_plan_year_matrix) * (average_claims .*
					plan_demands ./ plan_risk_scores));
			partial_firmclaims_partial_price = ac_partial_ownplans_partial_price .+ demand_partial_ac_partial_price;
			
		# Calculate marginal admin
		
			partial_firmadmin_partial_price = transpose(jacobian .* own_plan_year_matrix) * pmt[:,findall(pmt_fields .== "variable_admin")[1]];	
			
		# Calculate FOC Values
		
			partial_plan_demand_partial_price = diag(jacobian);
			
			foc_values = (transpose(foc_moment_indicator) * ((1 ./ partial_plan_demand_partial_price) .*
				(partial_firmrevenue_partial_price .- (1 .- pmt[:,findall(pmt_fields .== "reins_factor")[1]]) .* partial_firmclaims_partial_price .+ 
					(1 .- nora_flag) .* partial_firmtransfer_partial_price .- partial_firmadmin_partial_price) .* plan_demands)) ./
				(transpose(foc_moment_indicator) * plan_demands) .- foc_residuals;
		

			if pc_flag
				claims_reins = (1 .- pmt[:,findall(pmt_fields .== "reins_factor")[1]]) .* average_claims .* plan_demands;
				variable_admin_costs = pmt[:,findall(pmt_fields .== "variable_admin")[1]] .* plan_demands
				
				plan_rs_shares = plan_demands .* rs_cost_factor ./ sum(plan_demands .* rs_cost_factor);
				plan_ra_shares = plan_demands .* ra_cost_factor ./ sum(plan_demands .* ra_cost_factor);
				plan_transfers = (plan_rs_shares .- plan_ra_shares) .* annual_total_premiums;
				
				foc_values = transpose(foc_moment_indicator) * (plan_total_premiums .- claims_reins .+ (1 .- nora_flag) .* plan_transfers .- variable_admin_costs) ./
					(transpose(foc_moment_indicator) * plan_demands) .- average_plan_fixed_costs .- foc_residuals; 
			end

		return(foc_values)
	end	

	function ef_approach(base_premium,pmt,hdata,X_RS,X_CLMS,foc_residuals,mand_repeal_flag,voucher_flag,pc_flag,no_inertia_flag,no_uninsured_flag,ef_approach_flag)
	
		tolerance = 1;
		violations = ones(length(base_premium)) .* 100;
		new_premium = copy(base_premium);
		
		# Get all terms in Equilibrium Conditions
		while maximum(abs.(violations)) > tolerance
			println(cat(dims=1,plan_names[findmax(abs.(violations))[2]],violations[findmax(abs.(violations))[2]]))
			eq_results = update_foc_values(new_premium,pmt,hdata,X_RS,X_CLMS,foc_residuals,nora_flag,mand_repeal_flag,voucher_flag,pc_flag,no_inertia_flag,no_uninsured_flag,true);
			violations = eq_results[:,1];
			#CSV.write("ef_approach.csv",DataFrame(cat(dims=2,plan_year_names,new_premium,eq_results[:,1],eq_results[:,2])));
			new_premium = eq_results[:,2];
		end
		
		#CSV.write("ef_approach.csv",DataFrame(cat(dims=2,plan_year_names,base_premium,violations,new_premium)));
		
		return(new_premium)
	end
	
	function create_dhm_demands(pmt)
			
		dhm_demands = zeros(size(pmt)[1]);
		unique_rating_areas = unique(plans_pmt[!,:rating_area])
		for m in unique_rating_areas
			ra_indices = findall(plans_pmt[!,:rating_area] .== m);
			dhm_demands[ra_indices] .= .001/length(ra_indices) * sum(pmt[ra_indices,findall(pmt_fields .== "real_demand")[1]]);
		end
	
		return(dhm_demands)
	end
	
	# Calculate demand jacobian
	
	function get_voucher(hdata,X)
	
		# Create data frame
		silver_premiums = hdata[:,findall(hdata_fields .== "real_unsubs_premium")[1]];
		silver_premiums[nonsilver_indicators] = hdata[nonsilver_indicators,findall(hdata_fields .== "real_unsubs_premium")[1]] .+ 100000; # set high premium for everything but silver
		expdf = DataFrame(household_number = household_numbers,expvarall = exp.((X[uninsured_plans .== 0,:] * min_beta)/lambda),silver_premium = silver_premiums);
		
		# SLC Indices
		num_household_plans = combine(groupby(expdf,:household_number), :household_number => length)[!,:household_number_length];
		pop!(num_household_plans)
		num_household_plans = cat(dims=1,0,num_household_plans);
		starting_indices = cumsum(num_household_plans); # starting position for each household
		household_slc_indices = combine(groupby(expdf,:household_number), :silver_premium => determine_SLC_plan)[!,:silver_premium_determine_SLC_plan];
		slc_indices = starting_indices .+ household_slc_indices; # SLC indices in hdata object
		
		subsidy = max.(0,hdata[slc_indices,findall(hdata_fields .== "real_unsubs_premium")[1]] - households[!,:SLC_contribution]) .* 
					(households[!,:subsidized_members] .> 0); # each household's subsidy (0 if unsubsidized)
		
		return(subsidy[household_numbers])
	end
	